/*  sdcbrwpi.cpp  by K.Tsuru  */
// function ID 3553  DARDIX
/***********************************************************
pi by Borweins' method improved by by Takahashi and Kanada.
************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
SDouble BorweinsPi(){
    SDouble x(2.0);
    uint max_sz = x.MaxSize();
//It takes into irreducible maximum size mode.
//It confirms that your system has enough memory.

    SDouble a(x.Type(), max_sz), b(x.Type(), max_sz), y(x.Type(), max_sz), w(x.Type(), max_sz),
    d(x.Type(), max_sz);

    int c;
    d = Sqrt(2.0);
    a = 6 - 4*d;    y = 17 - 12*d;
    y.FixedPoint(1);
    int itrmax = 2*howpow2(max_sz);
    c = 0;
    do{
        d = ONE + RecQrrt(ONE-y);   // 1+(1-y)^(-1/4)
        y = ONE - 2/d;
        b = y*y;
        d = ONE+2*y+b;
        w = d*d;
        d = ONE+b;
        d = d*d;
        a = a*w-x*(w-d);
        if(y.IsMLT(ONE)) break; // y << 1.0 or not
        x *= 4;
        y = b*b; // This statement is lacked in Takahashi and Kanada's paper.
        c++;
    } while(c < itrmax);
    y.PointFree();
// free memory of figure[]
    b.SizeZero();   y.SizeZero();   w.SizeZero();   d.SizeZero();
    a = DReciprocal(a);
    if(c == itrmax) y.SetError(y.FATAL, "BorweinsPi", 3553);
    y.iterationCount = c;
    return a;
}
// obedient method
#if 0
SDouble BorweinsPi(){
    SDouble a, a1, y2, y4, y, d, pi, p(8.0);
    int c;
    d = Sqrt(2.0);
    a = 6 - 4*d;    y = d - ONE;
    y.FixedPoint(0);
    int itrmax = 2*howpow2(a.MaxSize());
    c = 0;
    do{
        y2 =y*y;    y4 = y2*y2; // y4 = y^4
        y4 = Qrrt(ONE - y4);    // (1-y4)^(1/4)
        y = (ONE-y4)/(ONE+y4);
        y2 = y*y;
        a1 = a*Dpow(ONE+y, 4) - p*(y*(ONE + y2)+y2);
        d = a1-a;
        if(d.IsMLT(a)) break;
        a = a1;
        p *= 4;
        c++;
    } while(c < itrmax);
    pi = DReciprocal(a);
    y.iterationCount = c;
    return pi;
}
#endif